#pragma once
#ifndef MATH_VECTOR_N_H_
#define MATH_VECTOR_N_H_

#include "../Common/GlobalDefines.h"
#include "../Common/ZException.h"
#include <cassert> // for throw divided by 0 exception
#include <algorithm>

GLOBAL_NAMESPACE_BEGIN

namespace MATH
{
    enum CVectorNType
    {
        CVN_Default,
        CVN_ZERO,
        CVN_MAX,
        CVN_MIN,
        CVN_RAND
    };

    template <typename T, int N>
    class CVectorN
    {
    public:
        typedef T (*DistanceFunc) (const CVectorN& v1, const CVectorN& v2);

    public:
        CVectorN(CVectorNType type = CVN_Default, DistanceFunc distFunc = defaultDistanceFunc);
        CVectorN(const CVectorN& rhs);
        CVectorN(double* v);
        CVectorN(T x, T y);
        CVectorN(T x, T y, T z);
        CVectorN(T x, T y, T z, T w);
        ~CVectorN(){}

        T           dotProduct( const CVectorN& rhs ) const;
        CVectorN    crossProduct(const CVectorN& rhs) const;
        CVectorN    operator+ ( CVectorN& rhs )   const;
        CVectorN    operator+ ( const CVectorN& rhs ) const;
        CVectorN    operator- ( CVectorN& rhs )   const;
        CVectorN    operator- ( const CVectorN& rhs ) const;
        CVectorN    operator* ( T k ) const;
        CVectorN    operator/ ( T t ) const;
        CVectorN&   operator= ( const CVectorN& rhs );
        CVectorN&   operator+=( const CVectorN& rhs );
        CVectorN&   operator-=( const CVectorN& rhs );
        CVectorN&   operator*=( T k );
        CVectorN&   operator/=( T t );
        T           operator[]( int ind ) const;
        T&          operator[]( int ind );

        T   x() const   { return m[0]; }
        T&  x()         { return m[0]; }
        T   y() const   { if (N >= 1) return m[1]; else ZFW_EXCEPTION(Exception::ERR_NOT_IMPLEMENTED, "Not support for N>=1.", "ZVectorN<T,N>::y()"); }
        T&  y()         { if (N >= 1) return m[1]; else ZFW_EXCEPTION(Exception::ERR_NOT_IMPLEMENTED, "Not support for N>=1.", "ZVectorN<T,N>::y()"); }
        T   z() const   { if (N >= 2) return m[2]; else ZFW_EXCEPTION(Exception::ERR_NOT_IMPLEMENTED, "Not support for N>=2.", "ZVectorN<T,N>::z()"); }
        T&  z()         { if (N >= 2) return m[2]; else ZFW_EXCEPTION(Exception::ERR_NOT_IMPLEMENTED, "Not support for N>=2.", "ZVectorN<T,N>::z()"); }

        void    unify();
        double  length() const;
        double  norm() const;
        double  norm2() const;

        operator T*()       { return m; }
        operator const T*() { return m; }

        T minimalElem() const;
        T maximalElem() const;

        static T defaultDistanceFunc(const CVectorN& v1, const CVectorN& v2);

    private:
        void init(CVectorNType type);

    private:
        T m[N];
        DistanceFunc dFunc;
    };

    // IO
    template <typename U, int N>
    std::ostream& operator<< (std::ostream& out, const CVectorN<U, N>& v)
    {
        for (int i = 0; i < N; i++)
        {
            out << v[i] << " ";
        }
        return out;
    }

    template <typename U, int N>
    std::istream& operator>> (std::istream& in, const CVectorN<U, N>& v)
    {
        for (int i = 0; i < N; i++)
        {
            in >> m[i];
        }
        return in;
    }

    // static members
    template <typename T, int N>
    T CVectorN<T, N>::defaultDistanceFunc(const CVectorN<T, N>& v1, const CVectorN<T, N>& v2)
    {
        return (v1-v2).norm2();
    }

    // implementation
    template <typename T, int N>
    CVectorN<T, N>::CVectorN(CVectorNType type, DistanceFunc func)
    {
        if (func==NULL)
            dFunc = defaultDistanceFunc;
        else
            dFunc = func;

        init(type);
    }

    template <typename T, int N>
    CVectorN<T, N>::CVectorN(const CVectorN<T, N>& rhs)
    {
        dFunc = rhs.dFunc;
        for (int i = 0; i < N; i++)
        {
            m[i] = rhs[i];
        }
    }

    template <typename T, int N>
    CVectorN<T, N>::CVectorN(double* v)
    {
        dFunc = defaultDistanceFunc;
        for (int i = 0; i < N; i++)
        {
            m[i] = v[i];
        }
    }

    template <typename T, int N>
    CVectorN<T, N>::CVectorN(T x, T y)
    {
        dFunc = defaultDistanceFunc;
        m[0] = x;
        if (N > 1) m[1] = y;
        for (int i = 2; i < N; i++)
        {
            m[i] = 0;
        }
    }

    template <typename T, int N>
    CVectorN<T, N>::CVectorN(T x, T y, T z)
    {
        dFunc = defaultDistanceFunc;
        m[0] = x;
        if (N > 1) m[1] = y;
        if (N > 2) m[2] = z;
        for (int i = 3; i < N; i++)
        {
            m[i] = 0;
        }
    }

    template <typename T, int N>
    CVectorN<T, N>::CVectorN(T x, T y, T z, T w)
    {
        dFunc = defaultDistanceFunc;
        m[0] = x;
        if (N > 1) m[1] = y;
        if (N > 2) m[2] = z;
        if (N > 3) m[3] = w;
        for (int i = 4; i < N; i++)
        {
            m[i] = 0;
        }
    }

    template <typename T, int N>
    CVectorN<T, N>& CVectorN<T, N>::operator=(const CVectorN<T, N>& rhs)
    {
        dFunc = rhs.dFunc;
        for (int i = 0; i < N; i++)
        {
            m[i] = rhs[i];
        }
        return (*this);
    }

    template <typename T, int N>
    void CVectorN<T, N>::init(CVectorNType type)
    {
        switch (type)
        {
        case CVN_MAX:
            for (int i = 0; i < N; i++)
            {
                m[i] = T(DBL_MAX);
            }
            break;
        case CVN_MIN:
            for (int i = 0; i < N; i++)
            {
                m[i] = T(-DBL_MAX);
            }
            break;
        case CVN_RAND:
            for (int i = 0; i < N; i++)
            {
                m[i] = T((T)rand()/RAND_MAX);
            }
            break;
        case CVN_Default:
        case CVN_ZERO:
        default:
            for (int i = 0; i < N; i++)
            {
                m[i] = T(0);
            }
            break;
        }
    }

    template <typename T, int N>
    T CVectorN<T, N>::dotProduct( const CVectorN<T,N>& v ) const
    {
        T ret = T(0);
        for (int i = 0; i < N; i++)
        {
            ret += m[i] * v[i];
        }
        return ret;
    }

    template <typename T, int N>
    double CVectorN<T, N>::norm2() const
    {
        double ret = 0;
        for (int i = 0; i < N; i++)
        {
            ret += m[i] * m[i];
        }
        return ret;
    }

    template <typename T, int N>
    double CVectorN<T, N>::norm() const
    {
        return sqrt(norm2());
    }

    template <typename T, int N>
    double CVectorN<T, N>::length() const
    {
        return norm();
    }

    template <typename T, int N>
    void CVectorN<T, N>::unify()
    {
        T len = length();
        if (len<DBL_EPSILON)
        {
            // this could be a random vector
            for (int i=0; i<N; i++)
            {
                m[i] = (T)(rand()%1000+20);
            }
            len=length();
            for (int i=0; i<N; i++)
            {
                m[i] /= len;
            }
        }
        else
        {
            for (int i = 0; i < N; i++)
            {
                m[i] /= len;
            }
        }
    }

    //template <typename T, int N>


    template <typename T, int N>
    CVectorN<T, N> CVectorN<T, N>::operator+(const CVectorN<T, N>& rhs) const
    {
        CVectorN ret;
        for (int i = 0; i < N; i++)
        {
            ret[i] = m[i] + rhs[i];
        }
        return ret;
    }

    template <typename T, int N>
    CVectorN<T, N> CVectorN<T, N>::operator+(CVectorN<T, N>& rhs) const
    {
        CVectorN ret;
        for (int i = 0; i < N; i++)
        {
            ret[i] = m[i] + rhs[i];
        }
        return ret;
    }

    template <typename T, int N>
    CVectorN<T, N> CVectorN<T, N>::operator-(const CVectorN<T, N>& rhs) const
    {
        CVectorN ret;
        for (int i = 0; i < N; i++)
        {
            ret[i] = m[i] - rhs[i];
        }
        return ret;
    }

    template <typename T, int N>
    CVectorN<T, N> CVectorN<T, N>::operator-(CVectorN<T, N>& rhs) const
    {
        CVectorN ret;
        for (int i = 0; i < N; i++)
        {
            ret[i] = m[i] - rhs[i];
        }
        return ret;
    }

    template <typename T, int N>
    CVectorN<T, N> CVectorN<T, N>::operator*(T k) const
    {
        CVectorN ret;
        for (int i = 0; i < N; i++)
        {
            ret[i] = m[i]*k;
        }
        return ret;
    }

    template <typename T, int N>
    CVectorN<T, N> CVectorN<T, N>::operator/(T t) const
    {
        assert(std::abs(t)>=DBL_EPSILON);
        CVectorN ret;
        for (int i = 0; i < N; i++)
        {
            ret[i] = m[i]/t;
        }
        return ret;
    }

    template <typename T, int N>
    CVectorN<T, N>& CVectorN<T, N>::operator+=(const CVectorN<T, N>& rhs)
    {
        for (int i = 0; i < N; i++)
        {
            m[i] += rhs[i];
        }
        return (*this);
    }

    template <typename T, int N>
    CVectorN<T, N>& CVectorN<T, N>::operator-=(const CVectorN<T, N>& rhs)
    {
        for (int i = 0; i < N; i++)
        {
            m[i] -= rhs[i];
        }
        return (*this);
    }

    template <typename T, int N>
    CVectorN<T, N>& CVectorN<T, N>::operator*=(T k)
    {
        for (int i = 0; i < N; i++)
        {
            m[i] *= k;
        }
        return (*this);
    }

    template <typename T, int N>
    CVectorN<T, N>& CVectorN<T, N>::operator/=(T k)
    {
        assert(std::abs(k)>=DBL_EPSILON);
        for (int i = 0; i < N; i++)
        {
            m[i] /= k;
        }
        return (*this);
    }

    template <typename T, int N>
    T& CVectorN<T, N>::operator[](int ind)
    {
        return m[ind];
    }

    template <typename T, int N>
    T CVectorN<T, N>::operator[](int ind) const
    {
        return m[ind];
    }

    template <typename T, int N>
    T CVectorN<T, N>::minimalElem() const
    {
        return (*std::min_element(m, m+N));
    }

    template <typename T, int N>
    T CVectorN<T, N>::maximalElem() const
    {
        return (*std::max_element(m, m+N));
    }

    template <typename T, int N>
    CVectorN<T, N> CVectorN<T, N>::crossProduct(const CVectorN<T, N>& rhs) const
    {
        if (N == 3)
        {
            CVectorN ret;
            ret[0] = m[1] * rhs[2] - m[2] * rhs[1];
            ret[1] = m[2] * rhs[0] - m[0] * rhs[2];
            ret[2] = m[0] * rhs[1] - m[1] * rhs[0];
            return ret;
        }
        else
        {
            ZFW_EXCEPTION(Exception::ERR_NOT_IMPLEMENTED, "Not support for cross producting vectors with dimension != 3. ", "CVectorN<T,N>::crossProduct");
        }
    }

    //// for class Color
    //// template specialization
    //template<> inline CVectorN<float, 4>::CVectorN(float x, float y, float z)
    //{
    //    m[0] = x;
    //    m[1] = y;
    //    m[2] = z;
    //    m[3] = 1.0;
    //    dFunc = NULL;
    //}

    typedef CVectorN<double, 2> vector2d;
    typedef CVectorN<double, 3> vector3d;
}

//typedef MATH::CVectorN<float, 4> Color;

GLOBAL_NAMESPACE_END

#endif//MATH_VECTOR_N_H_